Author

Victor L. Jardim

System and models

Real system

Two-Way Econometric Fixed Effects model

Group Mean Covariate Model # Richness models

Naive, mixed effects and fixed effects models

Code 
tar_load(alphamod)
alphamod <- alphamod %>% mutate(across(c(PC1_score, PC2_score, Depth, Fetch_max, Current_mean, T_mean, OM), scale)) %>% mutate(SY = paste(.$Site, .$Year, sep = "_"))


mod_naive <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean, data = alphamod)

mod_re <- lmer(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+(1|Site), data = alphamod)

mod_fe <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+Site+Year, data = alphamod)

mod_twfe <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+SY, data = alphamod)
Code 
export_summs(
  mod_naive,
  mod_re,
  mod_fe,
  mod_twfe,
  model.names = c("naive", "re", "fe", "twfe" ),
  scale = TRUE,
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "PC1" = "PC1_score",
    "PC2" = "PC2_score",
    # "PC1_site" = "PC1_site",
    # "PC2_site" = "PC2_site",
    # "PC1_sy" = "PC1_sy",
    # "PC2_sy" = "PC2_sy",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Fetch" = "Fetch_max",
    "Current" = "Current_mean",
    "Organic matter" = "OM",
    "Temperature" = "T_mean"),
  robust =  "HC1", error_pos = "right"
)
naiverefetwfe
PC12.46 *  (1.13)1.63   (1.33)-0.12   (1.15)0.20    (0.88)
PC24.26 ***(1.14)3.82 * (1.63)2.71   (1.45)2.08    (1.19)
PC1:PC20.43    (0.90)2.30 * (1.11)2.67 * (1.04)2.83 ***(0.82)
Fetch-2.52    (1.39)0.25   (3.27)13.90 **(4.39)13.33 ***(3.86)
Current-1.10    (1.16)-4.60   (3.45)-25.46 **(8.62)-18.04    (10.89)
Organic matter2.97 ** (0.95)4.07 **(1.23)2.30 * (0.98)3.62 ***(1.08)
Temperature2.44 ** (0.83)1.81 * (0.84)-1.81   (2.07)-16.84    (45.33)
N348           348          348          348           
DF339.00        337.00       319.00       222.00        
AIC2908.35        2873.04       2785.52       2712.69        
BIC2946.88        2915.42       2901.09       3201.92        
Marginal R2           0.26                            
Conditional R20.43        0.45       0.64       0.83        
Adjusted R20.42                  0.61       0.74        
F32.04                  20.55       8.94        
p-value0.00                  0.00       0.00        
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. Standard errors are heteroskedasticity robust. *** p < 0.001; ** p < 0.01; * p < 0.05.
Code 
(naivecoef <- coeftest(mod_naive, 
         vcov = vcov(mod_fe,
                       # cluster =  ~ Site,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive"))
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score              2.46       1.21     2.03  0.0432  naive
2 PC2_score              4.26       1.51     2.82  0.00513 naive
3 OM                     2.97       1.21     2.46  0.0145  naive
4 T_mean                 2.44       2.17     1.13  0.261   naive
5 Fetch_max             -2.52       4.80    -0.525 0.600   naive
6 Current_mean          -1.10       9.65    -0.114 0.910   naive
7 PC1_score:PC2_score    0.425      1.04     0.408 0.683   naive
Code 
(recoef <- tidy(mod_re) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive")) %>% 
  select(-effect:-group)
# A tibble: 7 × 5
  term                estimate std.error statistic model
  <chr>                  <dbl>     <dbl>     <dbl> <chr>
1 PC1_score              1.63      1.32     1.24   naive
2 PC2_score              3.82      1.57     2.44   naive
3 OM                     4.07      1.21     3.35   naive
4 T_mean                 1.81      0.837    2.17   naive
5 Fetch_max              0.249     3.08     0.0810 naive
6 Current_mean          -4.60      3.34    -1.38   naive
7 PC1_score:PC2_score    2.30      1.10     2.09   naive
Code 
coeftest(mod_fe, 
         vcov = vcovCL(mod_fe,
                       cluster =  ~ Site,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean"))
# A tibble: 7 × 5
  term                estimate std.error statistic  p.value
  <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
1 PC1_score             -0.118     1.10     -0.107 9.15e- 1
2 PC2_score              2.71      1.12      2.42  1.62e- 2
3 OM                     2.30      0.966     2.39  1.76e- 2
4 T_mean                -1.81      3.65     -0.496 6.20e- 1
5 Fetch_max             13.9       2.07      6.73  7.98e-11
6 Current_mean         -25.5       9.17     -2.78  5.83e- 3
7 PC1_score:PC2_score    2.67      1.09      2.44  1.51e- 2
Code 
(fecoef <- coeftest(mod_fe, 
         vcov = vcovCL(mod_fe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean")) |>
  mutate(model = "fe"))
# A tibble: 7 × 6
  term                estimate std.error statistic    p.value model
  <chr>                  <dbl>     <dbl>     <dbl>      <dbl> <chr>
1 PC1_score             -0.118     1.12     -0.106 0.916      fe   
2 PC2_score              2.71      1.02      2.66  0.00810    fe   
3 OM                     2.30      0.975     2.36  0.0187     fe   
4 T_mean                -1.81      3.06     -0.592 0.555      fe   
5 Fetch_max             13.9       2.94      4.73  0.00000334 fe   
6 Current_mean         -25.5      10.1      -2.52  0.0122     fe   
7 PC1_score:PC2_score    2.67      1.02      2.61  0.00944    fe
Code 
(twfecoef <- coeftest(mod_twfe, 
         vcov = vcovCL(mod_twfe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "twfe"))
# A tibble: 7 × 6
  term                estimate std.error statistic   p.value model
  <chr>                  <dbl>     <dbl>     <dbl>     <dbl> <chr>
1 PC1_score              0.202     1.29      0.156 0.876     twfe 
2 PC2_score              2.08      1.49      1.40  0.163     twfe 
3 OM                     3.62      0.824     4.39  0.0000177 twfe 
4 T_mean               -16.8      67.8      -0.248 0.804     twfe 
5 Fetch_max             13.3       3.23      4.13  0.0000522 twfe 
6 Current_mean         -18.0       8.96     -2.01  0.0453    twfe 
7 PC1_score:PC2_score    2.83      1.27      2.23  0.0269    twfe
Code 
coefs <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe")))

ggplot(coefs, aes(x = estimate, y = term, colour = model, shape = model))+ 
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.5, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

Naive, Grouped Mean Covariate, and Two-way grouped mean covariate models

Code 
alphamod <- alphamod %>%
  group_by(Site) %>% 
  mutate(PC1_site = mean(PC1_score),
         PC2_site = mean(PC2_score)
         ) %>% 
  ungroup() %>% 
  group_by(Year) %>% 
  mutate(PC1_y = mean(PC1_score),
         PC2_y = mean(PC2_score)) %>% 
  ungroup() %>% 
  group_by(Point) %>% 
  mutate(PC1_point = mean(PC1_score),
         PC2_point = mean(PC2_score)) %>% 
  ungroup()
  
mod_gmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + (1|Year), data = alphamod)
mod_fegmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + Year, data = alphamod)
mod_twgmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + PC1_y + PC2_y + Depth+OM+T_mean+Fetch_max + Current_mean + (1|Site) + (1|Year), data = alphamod)

export_summs(
  mod_naive,
  mod_gmcov,
  mod_fegmcov,
  mod_twgmcov,
  model.names = c("naive", "gmcov", "fegmcov", "twgmcov" ),
  scale = TRUE,
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "PC1" = "PC1_score",
    "PC2" = "PC2_score",
    "PC1_site" = "PC1_site",
    "PC2_site" = "PC2_site",
    "PC1_y" = "PC1_y",
    "PC2_y" = "PC2_y",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Fetch" = "Fetch_max",
    "Current" = "Current_mean",
    "Organic matter" = "OM",
    "Temperature" = "T_mean")
)
naivegmcovfegmcovtwgmcov
PC12.46    0.43   0.41   0.42   
(1.27)   (1.21)  (1.21)  (1.21)  
PC24.26 ***2.43   2.44   2.44   
(1.16)   (1.52)  (1.51)  (1.52)  
PC1_site       13.14 * 12.77 * 13.16 * 
       (5.65)  (5.78)  (5.67)  
PC2_site       11.63 * 11.69 * 11.62 * 
       (4.81)  (4.93)  (4.82)  
PC1_y                   2.70   
                   (3.53)  
PC2_y                   -0.88   
                   (3.64)  
PC1:PC20.43    3.02 **3.02 **3.02 **
(0.98)   (1.00)  (1.00)  (1.00)  
Fetch-2.52    10.60 * 10.82 * 10.63 * 
(1.43)   (4.43)  (4.46)  (4.43)  
Current-1.10    -4.76   -5.65   -4.87   
(1.15)   (4.58)  (4.74)  (4.61)  
Organic matter2.97 ** 3.30 **3.16 **3.19 **
(1.11)   (1.16)  (1.17)  (1.17)  
Temperature2.44 ** 0.01   -1.75   -0.10   
(0.85)   (1.72)  (2.17)  (1.81)  
N348       348      348      348      
DF339.00    334.00   324.00   332.00   
AIC2908.35    2799.41   2739.44   2794.70   
BIC2946.88    2853.34   2831.90   2856.33   
Marginal R2       0.36   0.46   0.37   
Conditional R20.43    0.70   0.71   0.71   
Adjusted R20.42                      
F32.04                      
p-value0.00                      
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05.
Code 
gmcoef <- tidy(mod_gmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "gmcov") %>% 
  select(-effect:-group)

fegmcoef <- tidy(mod_fegmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "fegmcov") %>% 
  select(-effect:-group)

twgmcoef <- tidy(mod_twgmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "twgmcov") %>% 
  select(-effect:-group)


coefs <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef, gmcoef, fegmcoef, twgmcoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe", "gmcov", "fegmcov", "twgmcov")))
ggplot(coefs, aes(x = estimate, y = term, colour = model, shape = model))+ 
    geom_vline(xintercept = 0, linetype="dashed")+
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.5, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

Density models

Code 
dens_naive <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean, data = alphamod)

dens_re <- lmer(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+(1|Site), data = alphamod)

dens_fe <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+Site+Year, data = alphamod)

dens_twfe <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+SY, data = alphamod)

dens_gmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + (1|Year), data = alphamod)

dens_fegmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + Year, data = alphamod)

dens_twgmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + PC1_y + PC2_y + Depth+OM+T_mean+Fetch_max + Current_mean + (1|Site) + (1|Year), data = alphamod)

(naivecoef <- coeftest(dens_naive, 
         vcov = vcov(dens_fe,
                       # cluster =  ~ Site,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive"))
# A tibble: 7 × 6
  term                estimate std.error statistic    p.value model
  <chr>                  <dbl>     <dbl>     <dbl>      <dbl> <chr>
1 PC1_score             0.134     0.0493     2.71  0.00706    naive
2 PC2_score            -0.0521    0.0614    -0.849 0.397      naive
3 OM                    0.196     0.0490     4.01  0.0000760  naive
4 T_mean                0.136     0.0879     1.55  0.123      naive
5 Fetch_max             0.228     0.195      1.17  0.244      naive
6 Current_mean         -0.305     0.392     -0.778 0.437      naive
7 PC1_score:PC2_score  -0.200     0.0423    -4.74  0.00000314 naive
Code 
(recoef <- tidy(dens_re) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "re") %>% 
  select(-effect:-group))
# A tibble: 7 × 5
  term                estimate std.error statistic model
  <chr>                  <dbl>     <dbl>     <dbl> <chr>
1 PC1_score             0.0541    0.0527     1.03  re   
2 PC2_score             0.125     0.0652     1.92  re   
3 OM                    0.151     0.0485     3.11  re   
4 T_mean                0.0735    0.0325     2.26  re   
5 Fetch_max             0.148     0.166      0.896 re   
6 Current_mean         -0.436     0.224     -1.95  re   
7 PC1_score:PC2_score   0.0360    0.0437     0.822 re
Code 
gmcoef <- tidy(dens_gmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "gmcov") %>% 
  select(-effect:-group)

fegmcoef <- tidy(dens_fegmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "fegmcov") %>% 
  select(-effect:-group)

twgmcoef <- tidy(dens_twgmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "twgmcov") %>% 
  select(-effect:-group)

(fecoef <- coeftest(dens_fe, 
         vcov = vcovCL(dens_fe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean")) |>
  mutate(model = "fe"))
Warning in sqrt(diag(se)): NaNs produced
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score             0.0295    0.0460     0.640  0.523  fe   
2 PC2_score             0.134     0.0592     2.26   0.0243 fe   
3 OM                    0.0864    0.0648     1.33   0.184  fe   
4 T_mean               -0.0869    0.124     -0.703  0.483  fe   
5 Fetch_max             0.246     0.257      0.957  0.339  fe   
6 Current_mean         -1.06      0.486     -2.18   0.0300 fe   
7 PC1_score:PC2_score   0.0309    0.0366     0.846  0.398  fe
Code 
(twfecoef <- coeftest(dens_twfe, 
         vcov = vcovCL(dens_twfe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "twfe"))
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score             0.0526    0.0472     1.12   0.266  twfe 
2 PC2_score             0.0954    0.0709     1.35   0.179  twfe 
3 OM                    0.139     0.0650     2.14   0.0337 twfe 
4 T_mean               -2.25      2.49      -0.905  0.367  twfe 
5 Fetch_max             0.190     0.330      0.576  0.565  twfe 
6 Current_mean         -0.596     0.592     -1.01   0.315  twfe 
7 PC1_score:PC2_score   0.0428    0.0384     1.12   0.266  twfe
Code 
coefdens <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef, gmcoef, fegmcoef, twgmcoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe", "gmcov", "fegmcov", "twgmcov")))

ggplot(coefdens, aes(x = estimate, y = term, colour = model))+
  geom_vline(xintercept = 0, linetype="dashed")+
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.6, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

RDA

Code 
tar_load(envcomp)
tar_load(bcdens)
tar_load(pal)

envcomp <- envcomp %>% mutate(Year = as.factor(Year))

list.hp4 <- list(Depth = envcomp %>% select(Depth),
                 Hydrodynamics = envcomp %>% select(Current_mean, Fetch_max),
                 Granulometry = envcomp %>% select(OM),
                 Temperature = envcomp %>% select(T_mean),
                 Complexity = envcomp %>% select(PC1_score, PC2_score),
                 Site = envcomp %>%select(Site),
                 Year = envcomp %>% select(Year))


rdamaerl.hp4 <- rdacca.hp::rdacca.hp(bcdens, list.hp4, var.part = TRUE)
upset_vp_victor(rdamaerl.hp4, cutoff = 0.005, int.var = "Complexity", pal = pal[1:8], title.cex = 22, axis.cex = 20, pch.size = 6, col.width = .8, effect.cex = 6)
Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code 
tar_read(upsetmain)

Code 
tar_read(triplotrda)